Modeling multiple time scales during glass formation with phase field crystals 
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The dynamics of glass formation in monatomic and binary liquids are studied numerically using a 
microscopic field theory for the evolution of the time-averaged atomic number density A stochastic 
framework combining phase field crystal free energies and dynamic density functional theory is shown 
to successfully describe several aspects of glass formation over multiple time scales. Agreement with 
mode coupling theory is demonstrated for underdamped liquids at moderate supercoolings, and a 
rapidly growing dynamic correlation length is found to be associated with fragile behavior. 
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A unified theoretical framework within which the glass 
transition may be understood does not currently exist. 
The most significant theoretical advances have been con- 
centrated near the early stages of slowing, leaving the in- 
termediate and late stages relatively poorly understood. 
Mode coupling theory [1 (MCT) and molecular dynam- 
ics [2], for example, have provided insight into the initial 
regime of slowing above the so-called crossover temper- 
ature Tc, but are ineffective when applied to the slower 
regimes that occupy roughly ten orders of magnitude in 
time between Tc and the glass transition temperature Tg. 

Time- or ensemble- averaged dynamic density func- 
tional theories [3]-[6] (DDFTs) have been proposed as a 
more efficient means of describing slow dynamics below 
Tc, but several key issues remain unresolved: which of 
the proposed equations of motion are most appropriate, 
whether the details of the free energy significantly influ- 
ence dynamics, and whether the detailed predictions of 
MCT can be reproduced and eventually improved upon 
by such theories. Mean-field DFT functionals are known 
to typically produce multivalley free energy landscapes 
in which an exponential number of aperiodic solid states 
coexist below a certain T [7-9 . However, the nature of 
the transition by which a liquid evolves toward and be- 
tween these aperiodic solid states upon quenching is influ- 
enced heavily by the microscopic dynamics and thus, in 
DDFT, the equation of motion employed. Approximate 
analytic results [5] indicate that two DDFT equations of 
motion may describe a MCT-type glass transition, but 
numerical simulations have confirmed only stretched ex- 
ponential decay and super- Arrhenius slowing in related, 
non-DDFT models [4 . Here we provide direct numeri- 
cal solutions for a candidate DDFT that considers both 
inertia and damping, and utilizes the simplest DFT free 
energy, the phase field crystal (PFC) class [10 . 

The dimensionless Helmholtz potential of a two com- 
ponent PFC system can be written [TT] 
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In this notation i = Aoi B^Ui = ni{f, t)^ni is the scaled 
time averaged number density of i particles, fii is the 
species average number density, is related to the liquid 
bulk modulus, qi sets the equilibrium distance between 
particles of the same species, qab sets that between A 
and B particles, and Hi^ Ai, and A2 are constants 

(see Refs. [TO^ for further discussion of how these param- 
eters relate to material properties). The terms multiplied 
by Hi discourage < and are the distinguishing fea- 
ture of the Vacancy or VPFC model [11 . A hard > 
cutoff enforces the physical interpretation of as a num- 
ber density and in doing so produces a range of nonlinear 
responses. The resulting solutions take the form of inter- 
acting time- averaged density peaks, with local regions of 
:^ representing unoccupied, or vacancy, sites. In 
terms of glass formation, one can obtain qualitatively 
similar results whether Hi = or Hi ^0, but only the 
Hi ^ case is presented in the following. 

The simplest dynamics conserving may be written 
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where t is dimensionless time, Di ~ T, and r]i is a gaus- 
sian stochastic noise variable with (772(^1,^1)77^(7^2,^2)) = 
V • V(5(fi — r2)(^(ti — ^2). A second option is the over- 
damped equation of DDFT, 
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where n{f^t) is generally set to ni{f^t) and 
{^^{nM>^{r2M)) = V- V[n(r,t)(5(ri-r2)(5(ti-t2)]. 
A third equation reintroduces some of the faster 
dynamics by also including an inertial or wave-like term. 
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where and /3j are constants [12). 
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FIG. 1: (color online). Glass formation and crystallization in 
the monatomic VPFC model, (a) F(q* , t) at various T where 
corresponds to the first peak maximum in (q), (b) (q) 
at same T as in (a), offset vertically by 0.5n with n = 0, 1, 
(c) Arrhenius plot of r* from F(q*,t), Inset: stretching ex- 
ponent /3* from fit to F{q\t) = exp [-(t/r*)^*] , (d) TTT 
diagram: quenches from T = 1.6 at various T, points de- 
note where crystallization occurred, ua — 0.15, — —0.9, 
qA = 1, WA = 0, UA = 1, Ha = 1500, To = 1000, aA = 1, 
13a = 0.01, fs = fAB = 0, Ax = 1.0, At = 0.02, V = 128^ 



Previous PFC simulations indicate that Eq. ([2| sup- 
ports metastable glassy states but in general produces 
a discontinuous, nucleation driven liquid to glass transi- 
tion [9]. Recent analyses of Eqs. ([3| and Q suggest that 
both may recover the class of MCT equations for the liq- 
uid dynamic correlators that successfully describe a wide 
range of glass forming behaviors [5]. Here we numer- 
ically investigate, without approximation, Eq. Q with 
n(r, t) = 1, where the inclusion of stochastic noise implies 
a time-averaged rather than ensemble- averaged interpre- 
tation of DDFT [13 . Equilibrium liquid states at high 
Di were quenched by lowering the stochastic noise am- 
plitude T = ToDi at a rate T, and the freezing transition 
was analyzed for onset of vitrification or crystallization. 

We begin with results for monatomic systems, out- 
lined in Fig. [T] For T > 1.6 the structure and 
dynamics are those of a normal liquid. The mea- 
sured intermediate scattering functions {Fij{q^t) = 
{Sni{q,0)Sn^{q^t)) / Fij{q^O)) decay exponentially, the 
corresponding average relaxation times show an Arrhe- 
nius T dependence, and the structure factors are char- 
acteristic of an equilibrium liquid state. The function 
S^{q) quantifies the structural correlations of the lo- 
calized peaks in the density field. We define Sfj{q) = 
{5nf {q)5n^* {q')) ^ where Snf{r) is a binary map of the 
positions of the local number density peaks. 

Below T 1.6, the liquid begins to show signs of 
nonequilibrium behavior and the onset of glass forma- 



tion. F{q^t) becomes increasingly stretched and be- 
gins to exhibit a shoulder, the average relaxation time 
briefly begins to grow with a super- Arrhenius T de- 
pendence, and a split second peak emerges in S^{q). 
But signs of glass formation persist only to the freezing 
temperature of the crystal, Tf. Below this point crys- 
tallization interrupts the apparent glass transition un- 
less the liquid is rapidly quenched well below Tf. The 
time-temperature-transformation {TTT) diagram shown 
in Fig. [l]i demonstrates this behavior. The profile of 
the nose feature is typical for materials with relatively 
marginal glass forming ability, such as metallic glasses. 
Since long-lived glassy states are not supported in the re- 
gion 0.6 < T < T/, one cannot study a gradual dynamic 
transition from liquid to glass. This behavior is expected 
for simple monatomic systems. 

Thus we proceed to binary liquids and outline in Fig. [2] 
the qualitative behavior of one such model system for a 
range of dynamic conditions, from highly underdamped 
{ai//3i = 100) to highly overdamped (ofi/ft = 0.01). The 
chosen model contains equal number densities of A and 
B atoms {ua = tib = 0.075), and the equilibrium spac- 
ing of A atoms is 20% smaller than that of B atoms 
{Ra/Rb = Qb/qa = 0.8). Only the NN correlations are 
plotted in Figs. |2^-c, where N denotes the full density 
field Ua + ^b- When damping dominates, an effectively 
stable glass with dynamics resembling those of a strong 
glass former is generated. The dynamic correlators are 
generally best fit as a single exponential decay for all ac- 
cessible T, with increased stretching as T is lowered, but 
any plateaus are absent or ill-defined in the ai/Pi < 1 
data. The relaxation times exhibit a nearly Arrhenius 
T-dependence over the entire accessible T-range. 

At the opposite extreme, when inertia dominates, a 
transition with dynamics characteristic of fragile liquids 
is generated. The dynamic correlators show both stretch- 
ing and clear plateauing as T is lowered, and the di- 
vergence of the relaxation time is well fit by the Vogel- 
Fulcher form (r = A exp [B / {T - Tq)]) . This diver- 
gence becomes increasingly super- Arrhenius at higher T 
as ai/Pi grows. The underdamped transition at this level 
of detail qualitatively resembles that described by MCT. 

The fragility of the PFC liquid therefore appears to be 
strongly linked with the balance of inertial and damping 
terms in Eq. Q, ai/ j3i. The degree of fragility is in turn 
correlated with the nominal spatial extent of cooperative 
dynamic behavior, which is set by an inherent length 
scale associated with the inertial term. This term gen- 
erates wave modes which propagate over a fixed length 
scale in a crystal before being damped, and the result- 
ing dynamic correlation length follows ^^^'^^^^ on/ Pi 
p[2j . In a normal liquid these correlations are largely 
suppressed by the low density and weak structural cor- 
relations, so that ^^"^^^^ <C But with greater 
supercooling, as the system becomes increasingly dense 
and solid-like, the inertial correlations survive over length 
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FIG. 2: (color online). Binary VPFC results for various 
damping conditions. FNN{q* ,t) at various T shown for ai//3i 
of (a) 0.01, (b) 1, and (c) 100. (d) Arrhenius plot of the struc- 
tural relaxation times r*^ , (e) the same data (NN only) scaled 
to clarify strong vs. fragile behaviors, with Vogel-Fulcher fits 
shown as solid lines and power law fits as dashed lines. In- 
set: (tti/ft)T^iv vs. S^n{q*)^ demonstrating deviation from 
scaling when inertial effects become large, (f) S^j^{q) for 
ai/l3i 100. UA^fiB = 0.075, -0.9, qa = 1, = 0.8, 

w^ = 0,u^ = 1, H^ 1500, Tq 1000, qab 8/9, Al 100, 
As = 0, Ax = 1.0, At = 0.025, V = 64^ 128^ or 256^ 



scales which likely approach ^^^^^^^ . Roughly, ^^^^ 

^crystal _ ^.^^^ aS T ^ [14j. 

The effects of this growing dynamic length scale are 
therefore especially prominent in highly underdamped 
systems, where its properties can be observed and quan- 
tified through finite size effects in F{q^t). For exam- 
ple, when ai/j3i = 100 finite size effects become numer- 
ically insurmountable below T 2^ 0.4. Measurements in- 
dicate that the average two point liquid static correlation 
length ^^^^^^ grows slowly, approximately as 1/T, while 
the dynamic correlation length grows more rapidly, as 
^uquid _ _ 2|^)-i±o.35 (ggg Yig. This indicates 

that the supercooled liquid exhibits heterogeneous dy- 
namics driven by strong inertial effects. Similar links 
between slowing dynamics and growing dynamic correla- 
tion lengths have been widely discussed [2j El [T5HT8] . 

A correlation between fragility and the length scale 
for cooperativity is consistent with existing interpreta- 
tions of strong and fragile liquids [19]. We also find a 
relevant link to recent experiments on colloidal glasses 
which demonstrate a transition from strong to fragile be- 



havior as the elastic properties of the colloidal particles 
become increasingly stiff [20 . When overdamped, Eq. Q 
describes a very soft, visco-elastic solid, while elastic stiff- 
ness and fragility both increase as damping is reduced. 
This is because ai ^ Vs ^ V^, where Vs is a sound speed 
and E is the relevant elastic modulus. Greater elastic 
stiffness should therefore correspond to reduced effective 
damping and, we expect, increased fragility This 
agrees with the trend found in Ref. [20 . 

Figure [s] shows simulation images of N{r) for the 
o^i/Pi = 100 system averaged over various times at 
T = 1.225, 0.541, and 0.420. Caging is apparent at short 
times for all T, but the long time averages at low T retain 
more of their original structure as the peaks exhibit less 
translational freedom. It is important to note that time 
averages are shown at equal multiples of each liquid's re- 
laxation time, not at equal so that time scales remain 
normalized as T is varied. The continuous but rapid de- 
cline in translational freedom as T is lowered signals a 
smooth transition from liquid-like to activated dynam- 
ics. This is consistent with the postulated crossover at 
Tc, below which relaxations are expected to be limited 
by increasingly rare, heterogeneously correlated cage es- 
cape events. This transition coincides with the emer- 
gence of the plateau in Fij{q^ t) and the split second peak 
in Sfj^{q) below T ^ 0.6, as shown in Fig. [2j 

Figure [4] shows data for the ai/j3i = 100 system rele- 
vant for comparison with the predictions of MCT. The 
nonergodicity parameter, fij{q) (height of the plateau 
in Fij{q,t)), is plotted in Fig. [4^ for T = 0.420. It 
follows the normal MCT behavior in which fij decays 
while oscillating in phase with Sij{q). Some of the dy- 
namic scaling behaviors predicted by MCT are tested in 
Fig. (4)3. Present results indicate that the von Schweidler 
scaling for late /^-relaxations (initial decay after plateau) , 
F{q^t) = f — B{t/r)^, is obeyed reasonably well over 2- 
3 orders of magnitude in time. The measured exponent 
b :^ 0.45±0.15 is comparable to typical values. Fits to the 
MCT critical decay power law (initial decay to plateau), 
F{q, t) = f ^ At-"", give a 0.3 ± 0.1. 

The late a-relaxations predicted by MCT are generally 
well- approximated by a stretched exponential decay. Our 
data are fit quite well by this form, as shown in Fig. [4]3, 
but with a stretching exponent P that decreases with T 
from approximately 1 to 0.6. MCT also predicts that the 
initial divergences of the fast and slow relaxation times 
follow power laws, ^ (T - Tc)-^/^^^) and r^, - (T - 
Tc)~^, respectively, where 7 = l/(2a) + 1/(26). Fits to 
these forms are shown in Fig. |2^, and though the Vogel- 
Fulcher fits are superior, the power laws are reasonably 
accurate through the early super- Arrhenius growth. 

These results confirm that the DDFT equation of mo- 
tion with inertia does in fact describe a glass transition, 
and that when damping is weak this transition strongly 
resembles both the structural glass transition observed 
for fragile glass formers as well as that predicted by MCT. 
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FIG. 3: (color online). Time- averaged density evolution in 
the supercooled binary liquid with ai/^i — 100. N{r) at 
T = 1.225 (top), 0.541 (middle), and 0.420 (bottom) averaged 
over the indicated multiples of each system's t^a^. A sub- 
cubic section of each cell has been removed to reveal a portion 
of the inner simulation. nA{r) time averages are black-orange, 
TiB(f) time averages are black-white. Bottom Right: Log-plot 
of static and dynamic correlation lengths vs. 1/(T — To). The 
dashed red line has a slope of 1.0. 
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FIG. 4: (color online). The highly underdamped a^/ft = 100 
model and MCT. (a) fij{q) and Sfj{q) /2.7 at T = 0.420, (b) 
FNN{q,t/T^]si) with various MCT scaling functions. Dashed: 
von Schweidler power law. Dotted: stretched exponential, 
Dashed-Dotted: critical decay power law. 



Our results are consistent with a picture in which fragility 
is driven by a large dynamic correlation length, which in 
some cases can be associated with large elastic moduli. 
A direct test of this association could be performed using 
colloidal systems such as those of Ref. [20 . By varying 
the degree of confinement, one could compare the relative 
magnitudes and growth rates of any dynamic correlation 
length as behavior is varied from strong to fragile. 
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